import pandas as pd
import numpy as np
from scipy import stats
import statsmodels.api as sm
import os
import matplotlib.pyplot as plt
from sklearn import linear_model
from sklearn.preprocessing import StandardScaler
import statsmodels.formula.api as smf
import warnings
from sklearn.model_selection import cross_val_score
from datetime import datetime
from sklearn.cluster import MiniBatchKMeans, AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram
from sklearn.preprocessing import StandardScaler, OneHotEncoder, PolynomialFeatures
from sklearn import metrics
from datetime import datetime
from itertools import product
import warnings

Откроем и преобразуем данные¶

all_data = []
for dic, folders, files in os.walk('/Users/shrlq/Documents/python_yandex/course 6/week_3/agg_data_'):
    for file in files:
        print (file)
        data_month = pd.read_csv(str('/Users/shrlq/Documents/python_yandex/course 6/week_3/agg_data_/' + file), ';', index_col=['tpep_pickup_hour'], 
                                 parse_dates=['tpep_pickup_hour'])
        all_data.append(data_month)
all_data = pd.concat(all_data)
all_data.sort_values(by = 'tpep_pickup_hour', inplace=True)
all_data = all_data.drop(['Unnamed: 0'], axis=1)
all_data.head()

df_03.csv
df_02.csv
df_06.csv
df_05.csv
df_04.csv

Задание 1¶

Составьте из данных о поездках прямоугольную таблицу так, чтобы по строкам было время, а по столбцам идентификатор ячейки (возьмите только те, которые были отобраны на второй неделе).

# получим временные ряды по каждому региону
data_ = all_data.pivot_table(columns=['Reg_ID'], aggfunc=sum, values=['count'], index=all_data.index)
data_.columns = data_.columns.droplevel()

# исключим ряды, среднее которых за май меньше 5.
only_may = data_.iloc[2160:-720, :]
data_ = data_.iloc[:, np.where(only_may.mean(axis=0) >= 5)[0]]
data_.head()

Запишем функции для последующих преобразований¶

# функция для добавления синусоидных и косинусоидных признаков
def sincos_feat_gen(data):
    T_ = np.arange(0, data.shape[0])+1
    week_seas = []
    df = pd.DataFrame()
    for n in np.arange(1, 7):
        week_sin = np.sin(T_*2*np.pi*n/168)
        week_seas.append(week_sin)
        week_cos = np.cos(T_*2*np.pi*n/168)
        week_seas.append(week_cos)
        
    df_= pd.DataFrame(data = week_seas).T
    sincos_col_names = ['sincos_' + str(x) for x in np.arange(0, len(np.arange(1, 7))*2)]
    df_.columns = sincos_col_names
    return (df_)

# функция для добавления дискретного признака дня недели 
def wd_feat_gen(data, name='weekday'):   
    data_ = pd.DataFrame()
    data_[name] = [str(d.weekday()) for d in data.index]
    return(data_)

# функция для кодирования категорийных признаков в бинарные
def cat_feat_gen(data, columns):
    # добавим признак категориальные признаки дней недели:
    temp_data = data.loc[:, columns] 

    # с помощью OneHotEncoder закодируем признаки по дням неделям:
    encoder = OneHotEncoder(drop='first', categories='auto')
    data_oh = encoder.fit_transform(temp_data).toarray()
    data_oh = pd.DataFrame(data_oh)
    
    data_oh.index = data.index
    cat_names = ['cat_' + str(x) for x in np.arange(0, len(data_oh.columns))]
    data_oh.columns = cat_names
    return(data_oh)

# функция для создания полиномиальных признаков
def polynomy_feat_gen(data, n=2):
    PolyNomy = PolynomialFeatures(degree = n)
    data_pn = PolyNomy.fit_transform(data)
    data_pn = pd.DataFrame(data_pn)
    return(data_pn)

# функция нахождения признака T+i в периоде с лагом i
def count_i_lag(count, i):
    count_ = []
    for n in np.arange(count.shape[0]):
        if n+i >=0 and n+i < count.shape[0]:
            count_.append(count[n+i])
        else:
            count_.append(np.nan)
    return(count_)

# функция для стандартизации выбранных признаков
def standardize_only_num_feat (data_train, data_predict):
    SC = StandardScaler()
    model = SC.fit(data_train)
    data_pr = SC.transform(data_predict)
    data_pr = pd.DataFrame(data_pr)
    return (data_pr)

def time_s_decompose(data):
    with plt.rc_context():
        plt.rc("figure", figsize=(15,5))
        return (sm.tsa.seasonal_decompose(data).plot())

# функции регресии с помощью машинного обучения с добавлением регуляризатора L1, чтобы исключить линейно-зависимые признаки. 
def auxilary_reg_residuals(X, y, alpha=0.5):
    SR = linear_model.SGDRegressor(penalty='l1', max_iter=1000, tol=1e-3, random_state=123, alpha=alpha)
    model = SR.fit(X = X, y = y)

    scor = round(cross_val_score(SR, cv=5, X = X, y = y, 
                    scoring='r2', ).mean(), 3)

    print ('Среднее значение R2: ' + str(scor))
    print (X.columns[np.argwhere(model.coef_ > 0)])
    data_feat = X.columns[np.argwhere(model.coef_ > 0)]

    # посчитаем остатки регрессии, чтобы по ним потом обучить модель ARIMA
    predictions = SR.predict(X)
    data_ = pd.DataFrame(predictions, columns=['predictions'])
    data_['resid'] = data_['predictions'].values - y.values
    data_['resid']
    
    data_.index = X.index
    return([data_, data_feat])

# функция дифференцирования остатов регресии
def dif_seas_plot(data, col_current='resid', col_dif = 'resid_dif', period=None, period_prev=0):
    %pylab inline
    print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data[col_current][period_prev:])[1])

    data[col_dif] = data[col_current] - data[col_current].shift(period)
    plt.figure(figsize(15,5))
    sm.tsa.seasonal_decompose(data[col_dif][period+period_prev:]).plot()
    print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data[col_dif][period+period_prev:])[1])
    return (data)

# функция, которая строит все возможные комбинации гиперпараметров
def param_comb(p, q, P, Q):
    param = list(product(list(range(p+1)), list(range(q+1)), list(range(P+1)), list(range(Q+1))))
    print ('Количетво комбинаций параметров: ' + str(len(param)))
    return(param)

# функция, которая перебирает все комбинации гиперпараметров, считает AIC и выбирает лучшую комбинацию гиперпараметров
def SARIMAX_calc(data, parameters, s, period_days=None):
    results = []
    best_results =[]
    best_aic = float("inf")
    warnings.filterwarnings('ignore')
    
    if period_days is None:
        period_days = data.shape[0]
        
    parameters = np.array(parameters)[np.argsort([sum(a) for a in parameters])]

    for param in parameters:
        #try except нужен, потому что на некоторых наборах параметров модель не обучается
        try:
            model=sm.tsa.statespace.SARIMAX(data.iloc[data.shape[0] - period_days:, 0], 
                                            exog = data.iloc[data.shape[0] - period_days:, 1:],
                                    order=(param[0], d, param[1]), 
                                    seasonal_order=(param[2], D, param[3], s)).fit(disp=-1)
        #выводим параметры, на которых модель не обучается и переходим к следующему набору
        except ValueError:
            print('wrong parameters:', param) 
            continue
        #выводим параметры в случае возникновения коллинеарности 
        except LinAlgError:
            print('Singular matrix:', param)
            continue
        
        aic = model.aic
        #сохраняем лучшую модель, aic, параметры 
        if aic < best_aic:
            best_model = model 
            best_aic = aic 
            best_param = param 
            best_results.append([param, model.aic])
            print ([param, model.aic])
        results.append([param, model.aic]) 
    return(results)

def datetime_separate(data, column, new_names):
    data[new_names[0]] = list(map(lambda x: datetime.strftime(x, '%Y-%m-%d'), data[column]))
    data[new_names[1]] = list(map(lambda x: x.hour, data[column])) 
    data.drop([column], axis=1, inplace=True)
    return(data)

Подготовим таблицы с признаками и целевыми значениями, разобьем выборки по периодам¶

# Сгенериреум вспомогательные принаки, общие для любых географических зон
wd_data = wd_feat_gen(data_)  # добавляем признак дня недели
sincos_df = sincos_feat_gen(data_) # добавляем синусные и косинусные прищнаки
cat_df = cat_feat_gen(wd_data, ['weekday']) # преобразуем категорийные признаки в бинарные

pn_data = pd.concat([wd_data, sincos_df, cat_df], axis=1) # объединяем все признаки
pn_df = polynomy_feat_gen(pn_data, n=2) # добавляем полиномиальные признаки

pn_df.describe()

# Для каждого ряда и часа конца истории запишем в таблицу целевые значения T+i, где i от 1 до 6
reg_ids = data_.columns
step = 6

targ_data_list=[]
for reg in reg_ids:
    reg_df = pd.DataFrame()
    reg_df['tpep_pickup_hour'] = data_.index
    reg_df['Reg_ID'] = reg
    #-- reg_df['count'] = list(data_.loc[:, reg])
    
    targets=[]
    for i in list(range(step)):
        target_ = count_i_lag(data_.loc[:, reg], i=i)
        targets.append(target_)
    targ_df = pd.DataFrame(targets, index=['targ_' + str(x) for x in list(range(6))]).T
     
    targ_data_ = pd.concat([targ_df, reg_df], axis=1) #объеденим полученные данные для ряда
    targ_data_list.append(targ_data_) #объеденим полученные данные по всем рядам

targ_data = pd.concat(targ_data_list)
targ_data.head()

# сохраним результат в файл
targ_data.to_csv('/Users/shrlq/Documents/python_yandex/course 6/week_4/targ_data_w4.csv', ';', index=False)

# разобьем целевые значения на обучающую, тестовую и итоговую выборки - для дальнейшего расчета ошибки
may_start = datetime.strptime('2016-04-30 23:00:00', '%Y-%m-%d %H:%M:%S')
may_end = datetime.strptime('2016-05-31 17:00:00', '%Y-%m-%d %H:%M:%S')

june_start = datetime.strptime('2016-05-31 23:00:00', '%Y-%m-%d %H:%M:%S')
june_end = datetime.strptime('2016-06-30 17:00:00', '%Y-%m-%d %H:%M:%S')

train_cond = np.where(targ_data['tpep_pickup_hour'] <= may_start)[0]
y_train = targ_data.iloc[train_cond, :] 
print ('Кол-во строк в обуч. выборке: ' + str(y_train.shape[0]))

test_cond = np.where((targ_data['tpep_pickup_hour'] >= may_start) & (targ_data['tpep_pickup_hour'] <= may_end))[0]
y_test = targ_data.iloc[test_cond, :]  
print ('Кол-во строк в тест. выборке: ' + str(y_test.shape[0]))

val_cond = np.where((targ_data['tpep_pickup_hour'] >= june_start) & (targ_data['tpep_pickup_hour'] <= june_end))[0]
y_val = targ_data.iloc[val_cond, :] 
print ('Кол-во строк в фин. выборке: ' + str(y_val.shape[0]))

Кол-во строк в обуч. выборке: 220320
Кол-во строк в тест. выборке: 75378
Кол-во строк в фин. выборке: 72930

Сделаем кластеризацию временных рядов географических зон¶

Задание 2¶

Перед проведением кластеризации стандартизуйте столбцы (вычтите выборочное среднее и поделите на выборочную дисперсию).

# сделаем обучающую выборку для кластеризации
reg_cond_train = np.where(data_.index <= may_start)[0]
reg_data_train = data_.iloc[reg_cond_train, :]
reg_feat_train = pn_df.iloc[reg_cond_train, :]

# сделаем стандартизацию самих рядов (необходимо для последующей кластеризации) и соответствующих признаков
reg_data_train_sc = standardize_only_num_feat(reg_data_train, reg_data_train)
reg_feat_train_sc = standardize_only_num_feat(reg_feat_train, reg_feat_train)
reg_feat_train_sc.index = reg_data_train.index

# сделаем стандартизацию признаков за более длительный период на основе модели, построенной на обуч. выборке
pn_df_sc = standardize_only_num_feat(reg_feat_train, pn_df)
pn_df_sc.index = data_.index

# объеденим целевые значения географических зон и общие признаки за весь период
reg_comv_df = pd.concat([data_, pn_df_sc], axis=1)
reg_comv_df.describe()

Задание 3¶

Кластеризуйте географические зоны по значениям стандартизованных рядов. Подберите число кластеров так, чтобы оно было не слишком большим, но ряды внутри кластеров имели похожую форму.

# сделаем кластеризацию на данных по апрель (обучающей выборке)
# подберем такое минимальное количество кластеров, чтобы std было не большим
k = 4
cl_data_T = reg_data_train_sc.T

BMKMeans = MiniBatchKMeans(n_clusters=k, batch_size=500, random_state=1001)
BMKMeans.fit(cl_data_T)
clusters_ = BMKMeans.labels_
#найдем центры кластеров
centers = BMKMeans.cluster_centers_
plt.hist(clusters_, bins=k)
# добавим в таблицу названия кластеров
cl_data_T['cluster'] = clusters_
print (np.std(cl_data_T['cluster'].value_counts()))
print (cl_data_T['cluster'].value_counts())

# сделаем сопоставление кластеров и регионов
# сопоставим полученные ряды с кластерами
mapping = pd.DataFrame([data_.columns, clusters_], index=['Reg_ID', 'cl']).T
mapping.head()

2.0615528128088303
3    27
0    27
1    26
2    22
Name: cluster, dtype: int64

Задание 4¶

В каждом кластере выберите наиболее типичный ряд (например, это может быть ряд, соответствующий центру кластера)

# визуализируем центры кластеров
x = np.arange(0, 400)
for i in np.unique(clusters_):
    plt.figure(figsize=(20, 5))
    for y in [a[1][:-1] for a in cl_data_T[cl_data_T['cluster'] == i].iterrows()]:
        plt.plot(x, y[:400]) 
        plt.plot(x, centers[i][:400], color='red', linewidth=6)
        plt.title('Кластер №' + str(i))
    plt.show()

# временные ряды центров кластеров объединим с признаками
cl_cent_data_train = pd.DataFrame(centers, columns=np.array(reg_data_train.index))
cl_df_comb_train = pd.concat([cl_cent_data_train.T, reg_feat_train_sc], axis=1)
cl_df_comb_train.head()

# построим декомпозицию
for index in cl_cent_data_train.index:
    time_s_decompose(cl_cent_data_train.loc[index, :])
    plt.title('Кластер №:' + str(index))

Построим регрессии для каждого кластера¶

Одновременно с этим проведем отбор регрессионных признаков¶

# построим регрессию на основе данных по апрель и посчитаем остатки
# сохраним отобранные с помощью регуляризации признаки
data_resid=[]
selected_features=[]
for ind, val in enumerate(cl_cent_data_train.index):
    print ('Кластер №' + str(val))
    temp_df_comb = 
    .iloc[:, [ind] + list(range(k, len(cl_df_comb_train.columns)))]
    temp_df_comb.columns = ['count'] + ['x_' + str(n) for n, i in enumerate(temp_df_comb.columns[1:])]
    mod = auxilary_reg_residuals(X = temp_df_comb.iloc[:, 1:], 
                                 y = temp_df_comb.iloc[:, 0],
                                 alpha=.02)
    data_resid_ind = mod[0]
    data_resid.append(data_resid_ind) # таблица с остатками
    selected_feat = mod[1] # названия признаков
    selected_features.append(temp_df_comb.loc[:, ['count'] + list(selected_feat)]) # таблица с признаками

Кластер №0
Среднее значение R2: 0.806
Index(['x_32', 'x_86', 'x_100', 'x_114', 'x_116', 'x_126', 'x_140', 'x_149',
       'x_159', 'x_163', 'x_173', 'x_174'],
      dtype='object')
Кластер №1
Среднее значение R2: 0.238
Index(['x_1', 'x_20', 'x_32', 'x_49', 'x_66', 'x_80', 'x_81', 'x_86', 'x_100',
       'x_107', 'x_108', 'x_114', 'x_134', 'x_174', 'x_175'],
      dtype='object')
Кластер №2
Среднее значение R2: 0.73
Index(['x_18', 'x_20', 'x_37', 'x_67', 'x_86', 'x_96', 'x_100', 'x_121',
       'x_126', 'x_128', 'x_138', 'x_152', 'x_160', 'x_164', 'x_174', 'x_180',
       'x_188', 'x_207'],
      dtype='object')
Кластер №3
Среднее значение R2: 0.72
Index(['x_32', 'x_49', 'x_80', 'x_86', 'x_100', 'x_107', 'x_114', 'x_140',
       'x_149', 'x_159', 'x_163', 'x_169', 'x_174', 'x_175', 'x_187'],
      dtype='object')

Найдем оптимальные гиперпараметры для ARIMA (зная остатки регрессии)¶

# сделаем сезонное дифференцирование с периодом 24 для каждого кластера
data_dif1=[]
for ind, val in enumerate(cl_cent_data_train.index):
    temp_data_r = dif_seas_plot(data_resid[ind], col_current='resid', col_dif='resid_dif', period=24)
    plt.show()
    data_dif1.append(temp_data_r)

Populating the interactive namespace from numpy and matplotlib
Критерий Дики-Фуллера: p=0.000001
Критерий Дики-Фуллера: p=0.000000

<Figure size 1080x360 with 0 Axes>

Populating the interactive namespace from numpy and matplotlib
Критерий Дики-Фуллера: p=0.000005
Критерий Дики-Фуллера: p=0.000000

<Figure size 1080x360 with 0 Axes>

Populating the interactive namespace from numpy and matplotlib
Критерий Дики-Фуллера: p=0.000000
Критерий Дики-Фуллера: p=0.000000

<Figure size 1080x360 with 0 Axes>

Populating the interactive namespace from numpy and matplotlib
Критерий Дики-Фуллера: p=0.000018
Критерий Дики-Фуллера: p=0.000000

<Figure size 1080x360 with 0 Axes>

# сделаем еще раз обычное дифференцирование
indices = np.array(range(0, k))
data_dif2=[]
for ind in indices:
    temp_data_r = dif_seas_plot(data_dif1[ind], col_current='resid_dif', col_dif='resid_dif2', period=1, period_prev=24)
    plt.show()
    data_dif2.append(temp_data_r)

Populating the interactive namespace from numpy and matplotlib
Критерий Дики-Фуллера: p=0.000000
Критерий Дики-Фуллера: p=0.000000

<Figure size 1080x360 with 0 Axes>

Populating the interactive namespace from numpy and matplotlib
Критерий Дики-Фуллера: p=0.000000
Критерий Дики-Фуллера: p=0.000000

<Figure size 1080x360 with 0 Axes>

Populating the interactive namespace from numpy and matplotlib
Критерий Дики-Фуллера: p=0.000000
Критерий Дики-Фуллера: p=0.000000

<Figure size 1080x360 with 0 Axes>

Populating the interactive namespace from numpy and matplotlib
Критерий Дики-Фуллера: p=0.000000
Критерий Дики-Фуллера: p=0.000000

<Figure size 1080x360 with 0 Axes>

Выберем векторы начальных приближений¶

# построим автокоррелограмы для подбора начальных гиперпараметров
for ind, val in enumerate(cl_cent_data_train.index):
    print ('Кластер: №' +  str(val))
    ax = plt.subplot(212)
    sm.graphics.tsa.plot_acf(data_dif2[ind].resid_dif2[24+1:].values.squeeze(), lags=100, ax=ax)
    pylab.show()

    ax = plt.subplot(212)
    sm.graphics.tsa.plot_pacf(data_dif2[ind].resid_dif2[24+1:].values.squeeze(), lags=100, ax=ax)
    pylab.show()

Кластер: №0

Кластер: №1

Кластер: №2

Кластер: №3

Задание 5¶

Для выбранных географических зон подберите на исходных рядах оптимальную структуру моделей — набор регрессионных признаков и значения гиперпараметров p,d,q,P,D,Q

#Выберем начальные приближения
Q = [3, 2, 1, 5] 
q = [1, 3, 3, 4] 
P = [5, 2, 1, 0] 
p = [3, 2, 1, 0]

from itertools import product
parameters = []
for n in range(k):
    param = param_comb(p[n], q[n], P[n], Q[n])
    parameters.append(param)

Количетво комбинаций параметров: 192
Количетво комбинаций параметров: 108
Количетво комбинаций параметров: 32
Количетво комбинаций параметров: 30

Для каждого прогнозируемого кластера подготовим таблицу с признаками¶

# для каждого кластера сделаем таблицу с целевым значением центра этого кластера и 
# набором признаков, который был ранее отобран для соответствующего кластера
k=4
cl_data_for_SARIMAX = []
for ind, reg in enumerate(cl_cent_data_train.index):
    cl_df = cl_df_comb_train.iloc[:, [ind] + list(range(k, len(cl_df_comb_train.columns)))]
    cl_df.columns = ['count'] + ['x_' + str(n) for n, i in enumerate(cl_df.columns[1:])]
    
    cl_df = cl_df.loc[:, selected_features[ind].columns]
    cl_data_for_SARIMAX.append(cl_df)
    
print (cl_data_for_SARIMAX[0].shape)
print (cl_data_for_SARIMAX[1].shape)
print (cl_data_for_SARIMAX[2].shape)
print (cl_data_for_SARIMAX[3].shape)

(2160, 13)
(2160, 16)
(2160, 19)
(2160, 16)

Найдем оптимальные гиперпараметры модели, которые будем использовать для прогнозирования рядов¶

result_for_clusters=[]
min_result_for_clusters=[]
for n in list(range(k)):
    result = SARIMAX_calc(cl_data_for_SARIMAX[n], parameters=parameters[n], s=24)
    result_min = result[np.argmin([x[1] for x in result])]
    result_for_clusters.append([n, result])
    min_result_for_clusters.append([n] + result_min)
min_result_for_clusters

# запишем лучшие результаты в отдельную переменную значения гиперпараметров p,d,q,P,D,Q
best_results = [(0, 1, 5, 3), 
                (1, 2, 0, 1), 
                (1, 3, 1, 1), 
                (0, 4, 0, 5)] 

d = 1
D = 1

Перейдем к прогнозированию рядов географических зон¶

Для каждого прогнозируемого ряда подготовим таблицу с признаками¶

# для каждого ряда сделаем таблицу с целевым значением этого ряда и 
# набором признаков, который был ранее отобран для соответствующего кластера
# за все даты до конца истории июне
R=102
data_for_SARIMAX = []
for ind, reg in enumerate(data_.columns):
    reg_df = reg_comv_df.iloc[:, [ind] + list(range(R, len(reg_comv_df.columns)))]
    reg_df.columns = ['count'] + ['x_' + str(n) for n, i in enumerate(reg_df.columns[1:])]
    
    clust_ind = int(mapping.loc[mapping['Reg_ID'] == reg, 'cl'])
    
    reg_df = reg_df.loc[:, selected_features[clust_ind].columns]
    data_for_SARIMAX.append(reg_df)

Для каждого ряда построим прогнозы на май и на июнь¶

# для каждого ряда посчитаем модель на основе данных по апрель, с выбранными в соответствии с кластером гиперпараметрами
# с помощью релевантной модели построим прогноз с мая на 6 часов для каждого часа каждого ряда
predictions_may = []
predictions_jun = []
datelist_ = data_for_SARIMAX[0].index
warnings.filterwarnings('ignore')

for ind, reg in enumerate(reg_ids):
    train_data = data_for_SARIMAX[ind].iloc[:-(31*24+30*24), :] 
    test_data = data_for_SARIMAX[ind]
    skip = 25

    clust = int(mapping.loc[mapping['Reg_ID'] == reg, 'cl'])
    clust_ind = np.arange(k)[np.array(cl_cent_data_train.index==clust)][0]
    param = best_results[clust_ind]

    # сделаем прогнозы для каждого ряда 
    cl_regids = list(mapping.loc[mapping.cl==clust, 'Reg_ID'])
    regids_origin_ind = np.where(np.isin(reg_ids, reg_ids_ord))[0]
    
    try: 
        model_train=sm.tsa.statespace.SARIMAX(train_data.iloc[:, 0], 
                                    exog = train_data.iloc[:, 1:],
                                    order=(param[0], d, param[1]), 
                                    seasonal_order=(param[2], D, param[3], 24)).fit(disp=-1)
        print (model_train.aic)


        # проверим остатки на предмет несмещенности, стационарности и отсутствия автокоррелированности:
        print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(model_train.resid[skip:], 0)[1])
        print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(model_train.resid[skip:])[1])
        print ('Среднее по остаткам: ' + str(model_train.resid[skip:].mean()))

        plt.figure(figsize(10,10))
        plt.subplot(211)
        plt.plot(model_train.resid[skip:])
        plt.ylabel(u'Residuals')
        plt.title('Регион №' + str(reg))
        plt.show()
        
        if (stats.ttest_1samp(model_train.resid[skip:], 0)[1] < .9) & (model_train.resid[skip:].mean() > 5):
            coerc = model_train.resid[skip:].mean()
        else: 
            coerc = 0

        # построим прогноз на 6 часов вперед и отобразим его на графике:
        prediction_ = model_train.predict(start=2159, end=2159+6, dynamic = True, exog = test_data.iloc[2159:2159+6, 1:])
        prediction_.index = datelist_[2159:2159+6+1]

        plt.figure(figsize(15,5))
        plt.plot(train_data[['count']][-150:], color='r', alpha = .3)
        plt.plot(prediction_[-150:], color = 'b')

        plt.ylabel('Count')
        plt.title('Регион №' + str(reg))
        pylab.show()

        # построим модель для прогнозов на май и июнь с параметрами как в обучающей модели
        model_test = sm.tsa.statespace.SARIMAX(
        test_data.iloc[:, 0], exog = test_data.iloc[:, 1:], order=(param[0], d, param[1]),
            seasonal_order=(param[2], D, param[3], 24))
        model_test_res = model_test.filter(model_train.params)

        # построим прогнозы на май с плавающим концом истории и запишем прогнозы:
        x=2159 # индекс соответствующий 30 апреля 23 ч
        for y in np.arange(0, 31*24-6+1):
            prediction_may = model_test_res.predict(start=x+y+1, end=x+y+6, 
                                               exog = test_data.iloc[x:x+y+6, 1:], dynamic=True)
            prediction_may = prediction_may + coerc
            predictions_may.append([reg] + [datelist_[x+y]] + list(prediction_may))           
        
        # построим прогнозы на июнь с плавающим концом истории и запишем прогнозы:
        x=2903 # индекс соответствующий 31 мая 23 ч
        for y in np.arange(0, 30*24-6+1):
            prediction_jun = model_test_res.predict(start=x+y+1, end=x+y+6, 
                                               exog = test_data.iloc[x:x+y+6, 1:], dynamic=True)
            prediction_jun = prediction_jun + coerc
            predictions_jun.append([reg] + [datelist_[x+y]] + list(prediction_jun))           
        print('Success_reg:', reg)
    
    # при ошибке записываем прогноз, пропускаем прогноз этого ряда:
    except ValueError:
        print('ValueError_reg:', reg)
        continue
    except LinAlgError:
        print('LinAlgError_reg:', reg)
        continue

17377.789390717302
Критерий Стьюдента: p=0.915206
Критерий Дики-Фуллера: p=0.000000
Среднее по остаткам: -0.03230892009371647

Success_reg: 1075
20149.152224045585
Критерий Стьюдента: p=0.944367
Критерий Дики-Фуллера: p=0.000000
Среднее по остаткам: 0.04046061216547875

Success_reg: 1076
19486.07579618243
Критерий Стьюдента: p=0.782699
Критерий Дики-Фуллера: p=0.000000
Среднее по остаткам: 0.13687967459492636

Success_reg: 1077
17874.798955977098
Критерий Стьюдента: p=0.808631
Критерий Дики-Фуллера: p=0.000000
Среднее по остаткам: -0.08250935338944958

Success_reg: 1125
20457.3764712813
Критерий Стьюдента: p=0.386465
Критерий Дики-Фуллера: p=0.000000
Среднее по остаткам: -0.5404213350249102

Success_reg: 1126
22143.871595079596
Критерий Стьюдента: p=0.297724
Критерий Дики-Фуллера: p=0.000000
Среднее по остаткам: -0.9632968881014289

Success_reg: 1127
23170.067259333126
Критерий Стьюдента: p=0.962857
Критерий Дики-Фуллера: p=0.000000
Среднее по остаткам: -0.054791822717927825

Success_reg: 1128
23910.617490386656
Критерий Стьюдента: p=0.453298
Критерий Дики-Фуллера: p=0.000021
Среднее по остаткам: -1.0461931992706412

Success_reg: 1129
24648.185865104242
Критерий Стьюдента: p=0.840979
Критерий Дики-Фуллера: p=0.000000
Среднее по остаткам: 0.332340429446051

Success_reg: 1130
20862.589776970563
Критерий Стьюдента: p=0.705348
Критерий Дики-Фуллера: p=0.000000
Среднее по остаткам: -0.25905813723286714

# преобразуем прогноз в таблицу и сохраним его
predictions_test = pd.DataFrame(predictions_may, columns=np.array(['Reg_ID', 
                                                                   'tpep_pickup_hour'] + list(np.arange(1, 7))))
predictions_test.head()

predictions_test.to_csv('predictions_test.csv', index=False)

Задание 6¶

Для каждой из R географических зон настройте на данных до апреля 2016 включительно модель ARIMA с гиперпараметрами, соответствующими кластеру этой зоны. Для каждого конца истории от 2016.04.30 23:00 до 2016.05.31 17:00 постройте прогноз на 6 часов вперёд и посчитайте в ноутбуке ошибку прогноза. Прогноз построен выше. Ниже - расчет ошибки.

На данных за май расчитаем ошибку (out_of_sample error)¶

# переименуем колонки в таблице y_test с целевыми значениями
y_test.rename(columns={'targ_0': '1', 'targ_1': '2', 'targ_2': '3', 'targ_3': '4', 'targ_4': '5', 'targ_5': '6', }, inplace=True)

# удалим колонку с данными на конец истории  в каждой таблице
#predictions_test.drop(['count'], axis=1, inplace=True)
#y_test.drop(['count'], axis=1, inplace=True)

# сделаем плоскую таблицу
predictions_test_melt = predictions_test.melt(id_vars=['Reg_ID', 'tpep_pickup_hour'])
predictions_test_melt.rename(columns={'variable': 'step', 'value':'prediction'}, inplace=True)
y_test_melt = y_test.melt(id_vars=['Reg_ID', 'tpep_pickup_hour'])
y_test_melt.rename(columns={'variable': 'step', 'value':'target'}, inplace=True)

# объеденим таблицы
data_final_test = y_test_melt.merge(predictions_test_melt, on=['Reg_ID', 'tpep_pickup_hour', 'step'])

# посчитаем ошибку
data_final_test['delta'] = np.absolute(data_final_test.prediction - data_final_test.target)
MAE_test = data_final_test[['delta']].mean(axis=0)
MAE_test

delta    38.828818
dtype: float64

Задание 7¶

Для каждой из RR географических зон настройте на данных до мая 2016 включительно модель ARIMA с гиперпараметрами, соответствующими кластеру этой зоны. Для каждого конца истории от 2016.05.31 23:00 до 2016.06.30 17:00 постройте прогноз на 6 часов вперёд и запишите все прогнозы в файл в формате id,y, где столбец id состоит из склеенных через подчёркивание идентификатора географической зоны, даты конца истории, часа конца истории и номера отсчёта, на который делается предсказание (1-6); столбец y — ваш прогноз.

Прогноз построен выше. Ниже форматирование файла для выполнения задания.

Запишем прогнозы на июнь в файл¶

from datetime import datetime
# преобразуем прогноз на июнь в таблицу и сохраним его
predictions_jun = pd.DataFrame(predictions_jun, columns=np.array(['Reg_ID', 
                                                                   'tpep_pickup_hour'] + list(np.arange(1, 7))))
predictions_jun.to_csv('predictions_jun.csv', index=False)

# приведем таблицу к формату для отправки задания
predictions_jun_melt = predictions_jun.melt(id_vars=['Reg_ID', 'tpep_pickup_hour'])

id_ = predictions_jun_melt['Reg_ID'].astype(str)
date_ = predictions_jun_melt['tpep_pickup_hour'].apply(lambda x: datetime.strftime(x, '%Y-%m-%d')).astype(str)
hour_ = predictions_jun_melt['tpep_pickup_hour'].apply(lambda x: datetime.strftime(x, '%H')).astype(int).astype(str)
step_ = predictions_jun_melt['variable'].astype(str)
y_ = predictions_jun_melt['value']

predictions_jun_submit = pd.DataFrame()
predictions_jun_submit['id'] = id_ + '_' + date_ + '_' + hour_ + '_' + step_
predictions_jun_submit['y'] = y_

predictions_jun_submit.to_csv('predictions_jun_submit.csv', index=False)
predictions_jun_submit.head(20)

Задание 8¶

Загрузите полученный файл на kaggle: https://inclass.kaggle.com/c/yellowtaxi. Добавьте в ноутбук ссылку на сабмишн.

Я не знаю, как дать ссылку на сабмишн. Я в настоящее время на 116 позиции в Leaderboard. https://inclass.kaggle.com/c/yellowtaxi/leaderboard

	Reg_ID	count
tpep_pickup_hour
2016-02-01	1841	0
2016-02-01	269	0
2016-02-01	988	0
2016-02-01	961	0
2016-02-01	668	0

	0	1	2	3	4	5	6	7	8	9	...	200	201	202	203	204	205	206	207	208	209
count	3624.0	3624.000000	3624.000000	3.624000e+03	3.624000e+03	3.624000e+03	3.624000e+03	3.624000e+03	3.624000e+03	3624.000000	...	3624.000000	3624.0	3624.0	3624.0	3624.000000	3624.0	3624.0	3624.000000	3624.0	3624.000000
mean	1.0	2.960265	0.013964	-3.463115e-03	1.496175e-03	2.830902e-03	2.868940e-03	-2.563841e-03	2.385256e-03	0.001626	...	0.145695	0.0	0.0	0.0	0.139073	0.0	0.0	0.139073	0.0	0.139073
std	0.0	1.996566	0.706065	7.081962e-01	7.066275e-01	7.077735e-01	7.071034e-01	7.072948e-01	7.074335e-01	0.706969	...	0.352849	0.0	0.0	0.0	0.346070	0.0	0.0	0.346070	0.0	0.346070
min	1.0	0.000000	-1.000000	-1.000000e+00	-1.000000e+00	-1.000000e+00	-1.000000e+00	-1.000000e+00	-9.972038e-01	-1.000000	...	0.000000	0.0	0.0	0.0	0.000000	0.0	0.0	0.000000	0.0	0.000000
25%	1.0	1.000000	-0.680173	-7.071068e-01	-6.801727e-01	-6.801727e-01	-7.071068e-01	-7.071068e-01	-6.801727e-01	-0.733052	...	0.000000	0.0	0.0	0.0	0.000000	0.0	0.0	0.000000	0.0	0.000000
50%	1.0	3.000000	0.037391	-1.960054e-15	1.042188e-15	1.284641e-15	1.530190e-15	-3.062236e-15	9.166282e-16	0.074730	...	0.000000	0.0	0.0	0.0	0.000000	0.0	0.0	0.000000	0.0	0.000000
75%	1.0	5.000000	0.707107	7.071068e-01	6.801727e-01	7.330519e-01	7.071068e-01	7.071068e-01	6.801727e-01	0.733052	...	0.000000	0.0	0.0	0.0	0.000000	0.0	0.0	0.000000	0.0	0.000000
max	1.0	6.000000	1.000000	1.000000e+00	1.000000e+00	1.000000e+00	1.000000e+00	1.000000e+00	9.972038e-01	1.000000	...	1.000000	0.0	0.0	0.0	1.000000	0.0	0.0	1.000000	0.0	1.000000

	targ_0	targ_1	targ_2	targ_3	targ_4	targ_5	tpep_pickup_hour	Reg_ID
0	19.0	6.0	5.0	4.0	4.0	12.0	2016-02-01 00:00:00	1075
1	6.0	5.0	4.0	4.0	12.0	28.0	2016-02-01 01:00:00	1075
2	5.0	4.0	4.0	12.0	28.0	32.0	2016-02-01 02:00:00	1075
3	4.0	4.0	12.0	28.0	32.0	91.0	2016-02-01 03:00:00	1075
4	4.0	12.0	28.0	32.0	91.0	54.0	2016-02-01 04:00:00	1075

	1075	1076	1077	1125	1126	1127	1128	1129	1130	1131	...	200	201	202	203	204	205	206	207	208	209
count	3624.000000	3624.000000	3624.000000	3624.000000	3624.000000	3624.000000	3624.000000	3624.000000	3624.000000	3624.000000	...	3624.000000	3624.0	3624.0	3624.0	3624.000000	3624.0	3624.0	3624.000000	3624.0	3624.000000
mean	65.663355	136.572572	86.867550	68.007450	186.457781	272.520419	300.571744	344.333609	411.111479	139.228753	...	0.003558	0.0	0.0	0.0	-0.015280	0.0	0.0	-0.015280	0.0	0.016884
std	40.929233	81.567132	58.227737	39.479084	97.408838	136.661118	169.044870	195.849721	240.270885	76.017506	...	1.003725	0.0	0.0	0.0	0.984441	0.0	0.0	0.984441	0.0	1.018049
min	0.000000	0.000000	0.000000	0.000000	0.000000	1.000000	1.000000	3.000000	1.000000	0.000000	...	-0.410891	0.0	0.0	0.0	-0.410891	0.0	0.0	-0.410891	0.0	-0.392232
25%	26.000000	54.750000	26.000000	34.000000	107.750000	170.750000	183.000000	217.000000	254.000000	86.000000	...	-0.410891	0.0	0.0	0.0	-0.410891	0.0	0.0	-0.410891	0.0	-0.392232
50%	70.000000	160.000000	95.000000	71.000000	199.000000	288.000000	293.500000	309.000000	359.000000	144.000000	...	-0.410891	0.0	0.0	0.0	-0.410891	0.0	0.0	-0.410891	0.0	-0.392232
75%	97.000000	202.000000	127.000000	94.000000	254.000000	361.000000	397.250000	446.000000	548.000000	188.000000	...	-0.410891	0.0	0.0	0.0	-0.410891	0.0	0.0	-0.410891	0.0	-0.392232
max	205.000000	369.000000	275.000000	247.000000	496.000000	681.000000	810.000000	1109.000000	1419.000000	476.000000	...	2.433737	0.0	0.0	0.0	2.433737	0.0	0.0	2.433737	0.0	2.549510

	0	1	2	3	1	2	3	4	5	...	200	204	207	209
2016-02-01 00:00:00	-0.832640	-0.476788	-0.115480	-1.299725	-1.493953	0.046399	1.431697	0.095356	1.418607	...	-0.410891	-0.410891	-0.410891	-0.392232
2016-02-01 01:00:00	-1.219843	-0.775130	-0.277184	-1.427938	-1.493953	0.099040	1.428722	0.200502	1.406771	...	-0.410891	-0.410891	-0.410891	-0.392232
2016-02-01 02:00:00	-1.508504	-1.029898	-0.539586	-1.492775	-1.493953	0.151533	1.423768	0.304469	1.387117	...	-0.410891	-0.410891	-0.410891	-0.392232
2016-02-01 03:00:00	-1.556172	-1.203369	-0.784256	-1.528458	-1.493953	0.203805	1.416842	0.406675	1.359756	...	-0.410891	-0.410891	-0.410891	-0.392232
2016-02-01 04:00:00	-1.537251	-0.982779	-0.754652	-1.433968	-1.493953	0.255783	1.407954	0.506549	1.324841	...	-0.410891	-0.410891	-0.410891	-0.392232

	Reg_ID	tpep_pickup_hour	1	2	3	4	5	6
0	1075	2016-04-30 23:00:00	69.653230	51.231108	37.588398	27.237721	30.120604	23.115830
1	1075	2016-05-01 00:00:00	51.913017	37.989502	27.503707	30.292545	23.228887	35.376788
2	1075	2016-05-01 01:00:00	30.438599	23.062213	27.347236	21.324957	34.124895	41.430555
3	1075	2016-05-01 02:00:00	12.713545	21.260077	17.288350	31.515520	39.714809	52.917380
4	1075	2016-05-01 03:00:00	21.911447	17.671491	31.769594	39.879049	53.025373	56.216273

	id	y
0	1075_2016-05-31_23_1	23.882802
1	1075_2016-06-01_0_1	18.235073
2	1075_2016-06-01_1_1	1.179826
3	1075_2016-06-01_2_1	0.751935
4	1075_2016-06-01_3_1	-0.173357
5	1075_2016-06-01_4_1	4.198240
6	1075_2016-06-01_5_1	18.692397
7	1075_2016-06-01_6_1	34.442428
8	1075_2016-06-01_7_1	63.430311
9	1075_2016-06-01_8_1	59.928386
10	1075_2016-06-01_9_1	51.957970
11	1075_2016-06-01_10_1	73.153557
12	1075_2016-06-01_11_1	100.909255
13	1075_2016-06-01_12_1	82.296624
14	1075_2016-06-01_13_1	118.053485
15	1075_2016-06-01_14_1	94.495766
16	1075_2016-06-01_15_1	80.568277
17	1075_2016-06-01_16_1	87.204845
18	1075_2016-06-01_17_1	97.259730
19	1075_2016-06-01_18_1	108.082505